%!TEX program = xelatex
%!TEX TS-program = xelatex
%!TEX encoding = UTF-8 Unicode

\documentclass[12pt,t,aspectratio=169,mathserif]{beamer}
%Other possible values are: 1610, 149, 54, 43 and 32. By default, it is to 128mm by 96mm(4:3).
%run XeLaTeX to compile.

\input{wang-slides-preamble.tex}

\begin{document}

\title{第8章：常微分方程}
%(1.1-1.2) 
%\institute{上海立信会计金融学院}
\author{JMS LQW}
%\date{2021年3月12日}

\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\begin{frame}[fragile=singleslide]{6.1.1. }
\begin{frame}{目录 }

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{enumerate}
\item[8.1.]  初值问题的标准形式
\item[8.2.]  基本思想，欧拉方法
\item[8.3.]  odeint 函数：范德波尔方程、洛伦兹方程
\item[8.4.]  线性边值问题，非线性边值问题：布拉图方程
\item[8.5.]  延迟微分方程：逻辑斯谛方程，麦克-格拉斯方程
\item[8.6.]  随机微分方程：Black-Scholes方程

\end{enumerate}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.1. 初值问题 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：将范徳波尔方程化为标准的一阶微分方程组，
$$\ddot{y} -\mu (1-y^2)\dot{y}+y=0, t\ge t_0, \,\,\, y(t_0)=y_0, \dot{y}(t_0)=v_0.$$

\vspace{0.5cm}

\item  答：先将未知函数与其一阶导函数写成一个向量，  
\begin{eqnarray*}
\vec{y}(t) = \begin{bmatrix} y(t) \\ \dot{y} (t) \end{bmatrix}, \hspace{0.2cm}
\vec{y}_0 = \begin{bmatrix} y_0 \\ v_0  \end{bmatrix}, \hspace{0.2cm}
f(\vec{y}) = \begin{bmatrix} \vec{y}[1] \\ \mu(1-\vec{y}[0]^2)\vec{y}[1]-\vec{y}[0] \end{bmatrix}. 
\end{eqnarray*}
则范德波尔方程的标准形式是 $$\dot{\vec{y}} (t) = f(\vec{y}(t),t), t\ge t_0, \,\,\, \vec{y}(t_0)=\vec{y}_0. $$

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.1. 范徳波尔方程 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  浏览布朗大学 Vladimir Dobrushkin 教授的 Mathematica 课程主页：

{\footnotesize \url{https://www.cfm.brown.edu/people/dobrush/am34/Mathematica/ch3/pol.html}
}

%\url{https://pi.math.cornell.edu/~templier/}

\item  问：简述 van der Pol 方程的研究历史。
 

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.2. 基本思想 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：设 $\lambda$ 是一个参数。考虑一阶微分方程的初值问题
 $$\dot{y} = \lambda(y-e^{-t}) - e^{-t},\,\, t\ge 0,\,\, y(0)=0.$$

\begin{enumerate}
\item  求该初值问题的解析解。
\item  对参数的不同取值，画出解析解的图像。
\item  用欧拉方法求数值解。
\end{enumerate}

\item  答：解析解为 $y=e^{-t} - e^{\lambda t}$. 



\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.2. }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{figure}
\centering
\includegraphics[height=0.7\textheight, width=0.7\textwidth]{pic/fig-8-2.png}
% \caption{ }
\end{figure}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.2. 前向欧拉方法 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：使用前向欧拉方法求解初值问题：$$\dot{y} = f(t,y) = \lambda(y-e^{-t}) - e^{-t},\,\, 0\le t\le 1,\,\, y(0)=0.$$

\item  答：自变量取步长为 $h=1/N$ 的等距网格 
$$(t_0,t_1,t_2,\cdots, t_n,\cdots, t_N) =(0, \frac{1}{N}, \frac{2}{N}, \cdots, \frac{n}{N}, \cdots, 1). $$ 
设数值解在这些格点的函数值为 $$(Y_0,Y_1,Y_2,\cdots, Y_n,\cdots, Y_N). $$
则前向欧拉方法的迭代公式为 $$Y_0=0,\,\, Y_{n+1} = Y_n + hf(t_n,Y_n). $$


\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.2. 单步误差 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：求单步误差。

\item  答：设解析解在所取网格上的函数值为 $$(y_0,y_1,y_2,\cdots,y_n,\cdots,y_N). $$
则单步误差为 $$\tau_n = \frac{y_{n+1}-y_n}{h} - f(t_n,y_n). $$

\item  问：什么是实际误差？

\item  答：实际误差是数值解与解析解的差，$$E_n = Y_n -y_n. $$

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.2. }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：实际误差与单步误差之间有什么样的递推关系？

\item  答：
\begin{eqnarray*}
E_{n+1} &=& Y_{n+1} - y_{n+1} \\ 
&=& Y_n +hf(t_n,Y_n) - ( h\tau_n + y_n + h f(t_n,y_n) ) \\
&=& E_n -  h\tau_n + h(f(t_n,Y_n)-f(t_n,y_n)) \\
&=& E_n -  h\tau_n + h\lambda(Y_n-y_n) \\
&=& (1+h\lambda) E_n -  h\tau_n. 
\end{eqnarray*}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.2. }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：%将$t_n$ 时刻的实际误差表示成 $t_0$ 时刻的实际误差，与单步误差的表达式。
将 $E_n$ 表示成 $E_0$ 与 $\tau_0,\tau_1,\cdots,\tau_{n-1}$ 的表达式。

\item  答：将所有迭代写出来，左右各自相加，抵消可得，
\begin{eqnarray*}
E_{n} &=& (1+h\lambda) E_{n-1} -  h\tau_{n-1} \\
(1+h\lambda)E_{n-1} &=& (1+h\lambda)^2 E_{n-2} -  h(1+h\lambda)\tau_{n-2} \\
(1+h\lambda)^2E_{n-2} &=& (1+h\lambda)^3 E_{n-3} -  h(1+h\lambda)^2\tau_{n-3} \\
\cdots && \cdots \\
(1+h\lambda)^{n-1}E_{1} &=& (1+h\lambda)^n E_{0} -  h(1+h\lambda)^{n-1}\tau_{0} \\
&\Rightarrow & \\
 E_{n} &=& (1+h\lambda)^n E_{0} - h \sum\limits_{m=1}^{n} (1+h\lambda)^m \tau_{m-1}. 
\end{eqnarray*}


\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.2. 刚性（stiffness） }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：什么是刚性？

\item  答：
\begin{enumerate}
\item  如果前向欧拉方法的数值解的误差是可控的，那么步长必须满足 $$|1+h\lambda|<1. $$
\item  当参数 $\lambda$ 很大时，需要设置很小的步长，这会导致数值计算非常慢。
\item  当这种现象发生时，称这个初值问题有刚性现象。
\end{enumerate}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.2. }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：分析前向欧拉方法的优缺点。

\item  答：
\begin{itemize}
\item  这是显式的数值方法。
\item  单步误差较小。
\item  可以局部改变步长，使得同时保持精度和稳定性。
\item  无法处理刚性问题。
\end{itemize}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.2. 后向欧拉方法 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：写出后向欧拉方法的迭代公式。

\item  答：在导数 $dy/dt=f(t,y)$ 的估计中，用后一时刻的未知的函数值代入，$$Y_{n+1} = Y_n + hf(t_n, Y_{n+1}). $$

\item  问：使用后向欧拉方法求解初值问题： 
\begin{eqnarray*}
\dot{y} = f(t,y) = \lambda(y-e^{-t}) - e^{-t},\,\, 0\le t\le 1,\,\, y(0)=0.
\end{eqnarray*}

\begin{enumerate}
\item  求实际误差的递推关系式。
\item  求实际误差的表达式。
\item  对给定参数 $\lambda$, 步长满足什么条件时，数值解是稳定的？
\end{enumerate}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.3. odeint 函数 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：对常微分方程初值问题，最著名的求解器是什么？

\begin{itemize}
\item  ODEPACK by Alan Hindmarsh at LLNL. 

\url{https://computing.llnl.gov/projects/odepack/software}

\item  scipy.integrate 模块中的 odeint 函数，背后是 LSODA 的 Fortran 程序。

\end{itemize}


\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.3.1. 理论背景 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：下述求解初值问题的 Python 函数的三个参数的含义分别是什么？
\begin{python}
y=odeint(func,y0,tvals)
\end{python}

\item  答：
\begin{itemize}
\item  第一个参数是标准形式的微分方程 $y'=f(y,t)$ 的右端的函数。
\item  第二个参数是 $t_0$ 时刻的初始值 $y_0$. 
\item  第三个参数是时间区间的一些分点 $[t_0,t_1,t_2,\cdots,t_n]$. 
\end{itemize}

\item  问：返回的结果是什么？

\item  答：返回初值问题的数值解 $[Y_0,Y_1,Y_2,\cdots,Y_n]$. 

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.3. }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：下述求解初值问题的 Python 函数的其余参数的含义分别是什么？
\begin{python}
y=odeint(func,y0,tvals,Dfun,atol,rtol)
\end{python}

\item  答：
\begin{itemize}
\item  第四个参数 Dfun 是函数 $f(y,t)$ 的雅可比矩阵。
\item  第五个参数 atol 是绝对误差 $\varepsilon_{abs}$ 的最小值，
%由局部误差小于绝对误差，
即 $|E_n|<\varepsilon_{abs}$.  %以此来调整步长 $h$. 
\item  第六个参数 rtol 相对误差 $\varepsilon_{rel}$ 的最小值，
%由局部误差小于相对误差，
即 $|E_n|<|Y_n|\varepsilon_{rel}$. %以此来调整步长。
\end{itemize}

\item  综合对绝对误差和相对误差的要求，选取步长 $h$ 的准则一般是 $$|E_n|<|Y_n|\varepsilon_{rel}+\varepsilon_{abs}. $$ 

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.3.2.2. 谐波振荡器 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：使用 odeint 求解初值问题 $$y'(t) = y(t), y(0)=1. $$

\begin{python}
import numpy as np
from scipy.integrate import odeint
f=lambda y,t:y
odeint(f,1,[0,1])
\end{python}

\item  结果输出未知函数 $y(t)$ 在 $t=0,1$ 的函数值，
\begin{python}
array([[1.        ],
       [2.71828193]])
\end{python}


\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.3.2.2. 简谐运动 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：使用 odeint 求解简谐运动问题 $$\ddot{y}(t) +\omega^2 y(t)=0, y(0)=1, \dot{y}=0. $$

\item  答：首先写成标准形式 
\begin{eqnarray*}
\vec{y} = \begin{bmatrix} y \\ \dot{y} \end{bmatrix}, \hspace{0.3cm} 
\dot{\vec{y}} = \begin{bmatrix} \dot{y} \\ \ddot{y} \end{bmatrix} = 
 \begin{bmatrix}  0&1 \\  -\omega^2 &0  \end{bmatrix}
 \begin{bmatrix} y \\ \dot{y} \end{bmatrix}, \hspace{0.3cm} 
\vec{y}_0 = \begin{bmatrix} 1 \\ 0 \end{bmatrix}. 
\end{eqnarray*}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.3.2.2. 简谐运动，代码 1/2 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{python}
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def rhs(Y,t,omega):
    y,ydot=Y
    return ydot,-omega**2*y

t_arr=np.linspace(0,2*np.pi,101)
y_init=[1,0]
omega=2.0
y_arr=odeint(rhs,y_init,t_arr,args=(omega,))
y,ydot=y_arr[:,0],y_arr[:,1]
\end{python}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.3.2.2. 简谐运动，代码 2/2 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{python}
fig=plt.figure()

ax1=fig.add_subplot(121)
ax1.plot(t_arr,y,t_arr,ydot)
ax1.set_xlabel('t')
ax1.set_ylabel('y and ydot')

ax2=fig.add_subplot(122)
ax2.plot(y,ydot)
ax2.set_xlabel('y')
ax2.set_ylabel('ydot')

plt.suptitle('Solution curve when omega =%2g' % omega)
fig.tight_layout()
fig.subplots_adjust(top=0.90)
\end{python}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.3.2.2. }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{figure}
\centering
\includegraphics[height=0.7\textheight, width=0.9\textwidth]{pic/fig-8-3-2-2.png}
% \caption{ }
\end{figure}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.3.2.3. 简谐运动问题的向量场 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：对简谐运动问题，画出相平面 $(y,\dot{y})$ 中的向量场，
$$\ddot{y}(t) +\omega^2 y(t)=0, y(0)=1, \dot{y}=0.$$

\item  代码 1/3: 
\begin{python}
import numpy as np
import matplotlib.pyplot as plt 
from scipy.integrate import odeint

def rhs(Y, t, omega):
    y,ydot=Y
    return ydot, -omega**2*y
\end{python}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.3.2.3. }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  代码 2/3: 

\begin{python}
t_arr=np.linspace(0, 2*np.pi, 101)
y_init=[1, 0]
omega=2.0

fig=plt.figure()
ax=fig.add_subplot(111)
y,ydot=np.mgrid[-3:3:21j, -6:6:21j]
u,v=rhs(np.array([y, ydot]), 0.0, omega)
mag=np.hypot(u, v)
mag[mag==0]=1.0
ax.quiver(y, ydot, u/mag, v/mag, color='blue')
\end{python}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.3.2.2. }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{figure}
\centering
\includegraphics[height=0.7\textheight, width=0.5\textwidth]{pic/fig-8-3-2-3.png}
% \caption{ }
\end{figure}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.3.2.2. 交互式画出轨线 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：点击相平面上的任意一点，画出从该点出发的轨线。

\item  代码 3/3: 

\begin{python}
# Enable drawing of arbitrary number of trajectories 
print('\n\n\nUse mouse to select each starting point')
print('Timeout after 30 seconds')
choice=[(0,0)]
while len(choice) > 0:
    y01 = np.array([choice[0][0], choice[0][1]]) 
    y = odeint(rhs, y01, t_arr, args=(omega,)) 
    plt.plot(y[:, 0], y[:, 1], lw=2)
    choice = plt.ginput()
print('Timed out!')
\end{python}


\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.3.3. 范德波尔振荡器 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：考虑范徳波尔方程，
$$\ddot{y} -\mu (1-y^2)\dot{y}+y=0, t\ge t_0, \,\,\, y((t_0)=y_0, \dot{y}(t_0)=v_0.$$

\begin{enumerate}
\item  写出标准形式 $\dot{\vec{y}} = f(\vec{y},t)$ 的右端函数的雅可比矩阵。
\item  使用 odeint 求解这个方程。
\end{enumerate} 

\item  答：标准形式的微分方程 
\begin{eqnarray*}
\begin{bmatrix} \dot{y} \\ \dot{z}  \end{bmatrix} = \begin{bmatrix} f(y,z) \\ g(y,z)  \end{bmatrix} 
\end{eqnarray*}
的右端函数的雅可比矩阵为 
\begin{eqnarray*}
\begin{bmatrix} f_y & f_z \\ g_y & g_z \end{bmatrix}. 
%\begin{bmatrix} \frac{\partial f}{\partial y} & \frac{\partial f}{\partial z} \\ \frac{\partial g}{\partial y} & \frac{\partial g}{\partial z} \end{bmatrix}.
\end{eqnarray*}


\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.3.3. }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  代码 1/2: 

\begin{python}
import numpy as np
from scipy.integrate import odeint

def rhs(z,t,mu):
    return [ z[1], mu*(1-z[0]**2)*z[1]-z[0] ]

def jac(z,t,mu):
    return[ [0,1], [-2*mu*z[0]*z[1]-1,mu*(1-z[0]**2)] ]
\end{python}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.3.3. }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  代码 2/2: 

\begin{python}
mu=1.0
t_final = 15.0 if mu<10 else 4.0*mu
n_points = 1001 if mu<10 else 1001*mu
t=np.linspace(0,t_final,n_points)
y0=np.array([0.1,0.0])
y,info=odeint(rhs,y0,t,args=(mu,),Dfun=jac,full_output=True)

print('mu=%g, number of Jacobian calls is %d' %
      (mu,info['nje'][-1]) )

import matplotlib.pyplot as plt
plt.plot(y[:,0],y[:,1])
\end{python}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.3.3. }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{figure}
\centering
\includegraphics[height=0.7\textheight, width=0.5\textwidth]{pic/fig-8-3-3.png}
% \caption{ }
\end{figure}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.3.4. 洛伦兹方程 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：洛伦兹方程是研究地球天气的数学模型简化而来的一个微分方程， 
\begin{eqnarray*}
\dot{x} &=& \sigma (y - x), \\
\dot{y} &=& \rho x - y - xz, \\
\dot{z} &=& xy - \beta z.
\end{eqnarray*}

设 $\sigma=10, \rho=29, \beta=8/3$. 求解这个微分方程。

\item  代码1/3: 

\begin{python}
import numpy as np
from scipy.integrate import odeint

def rhs(u,t,beta,rho,sigma):
    x,y,z=u
    return [sigma*(y-x), rho*x-y-x*z, x*y-beta*z ]
\end{python}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.3.4. }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  代码2/3: 

\begin{python}
sigma=10.0
beta=8.0/3.0
rho1=29.0
rho2=28.8

u01=[1.0,1.0,1.0]
u02=[1.0,1.0,1.0]

t=np.linspace(0.0, 50.0, 10001)
u1=odeint(rhs,u01,t,args=(beta,rho1,sigma))
u2=odeint(rhs,u02,t,args=(beta,rho2,sigma))

x1,y1,z1=u1[:,0],u1[:,1],u1[:,2]
x2,y2,z2=u2[:,0],u2[:,1],u2[:,2]
\end{python}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.3.4. }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  代码3/3: 

\begin{python}
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

fig=plt.figure()
ax=Axes3D(fig)
ax.plot(x1,y1,z1,'b-')
#ax.plot(x2,y2,z2,'g:')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.set_title('Lorenz equations with rho=%g '%rho1)
\end{python}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.3.4.  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{figure}
\centering
\includegraphics[height=0.7\textheight, width=0.7\textwidth]{pic/fig-8-3-4.png}
% \caption{ }
\end{figure}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.3.4. }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：求洛伦兹方程的平衡点。

\item  问：将上图用动画呈现出来。

\item  问：有什么关于洛伦兹方程的文献？

\begin{itemize}
\item  Sparrow, C. The Lorenz Equations, Springer, 1982.

\end{itemize}


\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.4. 两点边值问题 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：什么是两点边值问题？

\item  答：设 $y(x)$ 是定义在区间 $a\le x\le b$ 的未知函数，符合如下常微分方程和边界条件
$$\ddot{y} = f(x,y,\dot{y}), \,\,\, y(a)=A,y(b)=B. $$

\item  问：介绍一些常微分方程边值问题的参考文献？

\begin{enumerate}
\item  Ascher, U.M., Mattheij, R.M.M. and Russell, R.D. {\color{blue}Numerical Solution of Boundary Value Problems for Ordinary Differential Equations.} SIAM, 1995. 
\item  Ascher, U.M., Mattheij, R.M.M., Russell, R.D. and Petzold, L.R. {\color{blue}Computer Methods for Ordinary Differential Equations and Differential Algebraic Equations.} SIAM, 1998. 
\end{enumerate}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.4.1.  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：对下述线性微分方程的边值问题，解释打靶法的思路。
$$\ddot{y} = - h(x)y(x), \,\, 0\le x\le 1, \,\, y(0)=A,y(1)=B. $$

\item  答：
\begin{enumerate}

\item  首先求解两个初值问题：
\begin{eqnarray*}
\ddot{y}_1 = - h(x)y_1(x), \,\, y_1(0)=1,\dot{y}_1(0)=0. \\ 
\ddot{y}_2 = - h(x)y_2(x), \,\, y_2(0)=0,\dot{y}_2(0)=1. 
\end{eqnarray*}

\item  然后设边值问题的解为 $y(x)=C_1y_1(x) + C_2y_2(x)$,则边值条件可得
\begin{eqnarray*}
C_1y_1(0) + C_2y_2(0) = A, \\
C_1y_1(1) + C_2y_2(1) = B. 
\end{eqnarray*}

\end{enumerate}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.4.2. 边值问题的公式化 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：将下述特征值问题化为标准形式，$$u''(x) + \lambda u(x) =0, \,\, u(0)=u(1)=0, \,\, \int_0^1 u^2(s)ds=1.$$

\item  答：引入辅助变量 $$v(x)=\lambda, \,\, w(x) =\int_0^x u^2(s)ds,$$ 则微分方程和边界条件分别为
\begin{eqnarray*}
&& u''(x) = -u(x)v(x), \,\, v'(x)=0, \,\, w'(x) =u^2(x), \\
&& u(0)=0,\,\, w(0)=0,\,\, u(1)=0,\,\, w(1)=1. 
\end{eqnarray*}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.4.2.}
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：边值问题的数值求解的思路与文献？


\begin{enumerate}
\item  Bader, G. and Ascher, U.M. A new basis implementation for a mixed order boundary value ODE solver. SIAM J. Sci. Stat. Comp. 8, 483-500.  

\item  Bader 和 Ascher 的 Fortran 程序 colnew. %, scikits.bvp1lg. 

\url{http://pv.github.io/scikits.bvp1lg/colnew.html}

\end{enumerate}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.4.3. 两点边值问题的简单例子：简谐振动 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：使用 colnew 求解微分方程边值问题 $$\ddot{u}(x) + u(x) = 0, \,\, 0\le x\le \pi, \,\, u(0)=0, \dot{u}(\pi)=1. $$



\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.4.3. 另一个例子：贝塞尔方程 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：求解微分方程边值问题：
\begin{eqnarray*}
&& \ddot{u}(x) = - \frac{\dot{u}(x)}{x} + \left( \frac{\nu^2}{x^2}-1\right) u(x), \,\, 1\le x\le 10, \\
&& u(1)=J_\nu(1), \,\, u(10)=J_\nu(10), 
\end{eqnarray*}
其中 $J_\nu$ 是第一类贝塞尔函数 $$J_\nu(x) = \sum\limits_{k=0}^{\infty} \frac{(-1)^k}{\Gamma(\nu+k+1)\Gamma(k+1)}\left(\frac{x}{2}\right)^{2k+\nu}. $$

%\item  问：求解配对的微分方程边值问题：
%\begin{eqnarray*}
%&& \dot{v}(x) = c^{\nu+1}u(x), \,\, 1\le x\le 10, \\
%&& v(5) = 5^{\nu+1}J_{\nu+1}(5). 
%\end{eqnarray*}

\item  这个例子来自 bvp1lg 的官网：
\url{http://pv.github.io/scikits.bvp1lg/examples.html}

%\begin{python}
%scikits.bvp1lg.examples.bessel_colnew()
%\end{python}


\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.4.3.  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  Bessel 方程代码 1/4: 

\begin{python}
import numpy as np
import scipy.special as special
import scikits.bvp1lg.colnew as colnew

nu = 3.4123; degrees = [2, 1]
def fsub(x, z):
    u, du, v = z     # it's neat to name the variables
    return np.array([-du/x + (nu**2/x**2 - 1)*u, 
                     x**(nu+1) * u])
def dfsub(x, z):
    u, du, v = z
    zero = np.zeros(x.shape)
    return np.array([[(nu**2/x**2 - 1), -1/x, zero],
                      [       x**(nu+1), zero, zero]])
\end{python}


\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.4.3.  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  Bessel 方程代码 2/4: 

\begin{python}
boundary_points = [1, 5, 10]

def gsub(z):
    u, du, v = z
    return np.array([u[0] - special.jv(nu,   1),
            v[1] - 5**(nu+1) * special.jv(nu+1, 5),
            u[2] - special.jv(nu,   10)])

def dgsub(z):
    return np.array([[1, 0, 0],
                     [0, 0, 1],
                     [1, 0, 0]])
\end{python}


\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.4.3.  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  Bessel 方程代码 3/4: 

\begin{python}
tol = [1e-5, 0, 1e-5]
solution = colnew.solve(
        boundary_points, degrees, fsub, gsub,
        dfsub=dfsub, dgsub=dgsub,
        is_linear=True, tolerances=tol,
        vectorized=True, maximum_mesh_size=300)

solution.nmesh

x = np.linspace(1, 10, 101)
np.allclose(solution(x)[:,0], special.jv(nu, x),
            rtol=1e-4, atol=1e-8)
\end{python}


\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.4.3.  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  Bessel 方程代码 4/4: 

\begin{python}
np.allclose(solution(x)[:,2], 
            x**(nu+1)*special.jv(nu+1, x),
            rtol=1e-4, atol=1e-8)
            
import matplotlib.pyplot as plt
plt.plot(solution.mesh, 
         solution(solution.mesh)[:,2], '.',
         x, x**(nu+1)*special.jv(nu+1, x), '-')
         
#plt.show()
\end{python}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.4.3. 贝塞尔方程的解}
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{figure}
\centering
\includegraphics[height=0.7\textheight, width=0.7\textwidth]{pic/fig-8-4-3.png}
% \caption{ }
\end{figure}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.4.4. 线性特征值问题 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：求解特征值问题
\begin{eqnarray*}
\left\{ \begin{array}{rcl}
u''(x) &=& -u(x)v(x), \\ 
v'(x) &=& 0, \\ 
w'(x) &=& u^2(x),\\
\end{array}\right. 
\hspace{0.3cm}
\left\{ \begin{array}{rcl}
u(0) &=& 0, \\ 
w(0) &=& 0, \\ 
u(1) &=& 0, \\ 
w(1) &=& 1. \\
\end{array}\right. 
\end{eqnarray*}


\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.4.5. 非线性边值问题 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：使用 colnew 求解布拉图方程 $$\ddot{u}(x) + \lambda e^{u(x)} = 0, \,\, 0<x<1, \,\, \lambda > 0, \,\, u(0) = u(1) = 0. $$

\item  答：
\begin{enumerate}
\item  解析解： 设 $v(x) =\dot{u}(x)$. 将原方程转化成关于 $v$ 的微分方程，然后求解。
\item  数值解： 设 $v(x)=\lambda$, $w(x) = \int_0^x u^2(s)ds$, 将原方程写成另一种标准形式。

\end{enumerate}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.4.5. }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：使用 {\verb+scipy.integrate.solve_bvp+} 求解布拉图方程？

{\tiny 
\url{https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_bvp.html}
}

\item  答：记 $\vec{y}[0] =y$, $\vec{y}[1]=\dot{y}$, 原方程化成如下标准形式
\begin{eqnarray*}
\frac{d}{dx}\vec{y} = \vec{f}(x,\vec{y}) \,\, \Rightarrow \,\,
\frac{d}{dx} \begin{bmatrix} y \\ \dot{y}  \end{bmatrix} = \begin{bmatrix} \dot{y} \\ -\lambda e^{y}  \end{bmatrix}. 
\end{eqnarray*}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.4.5. }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item 代码 1/2:  

\begin{python}
import numpy as np

def fun(x, y):
    return np.vstack((y[1], -np.exp(y[0])))

def bc(ya, yb):
    return np.array([ya[0], yb[0]])

x = np.linspace(0, 1, 5)

y_a = np.zeros((2, x.size))
y_b = np.zeros((2, x.size))
y_b[0] = 3
\end{python}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.4.5. }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item 代码 2/2:  

\begin{python}
from scipy.integrate import solve_bvp
res_a = solve_bvp(fun, bc, x, y_a)
res_b = solve_bvp(fun, bc, x, y_b)

x_plot = np.linspace(0, 1, 100)
y_plot_a = res_a.sol(x_plot)[0]
y_plot_b = res_b.sol(x_plot)[0]

import matplotlib.pyplot as plt
plt.plot(x_plot, y_plot_a, label='y_a')
plt.plot(x_plot, y_plot_b, label='y_b')
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
\end{python}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.4.5. }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{figure}
\centering
\includegraphics[height=0.7\textheight, width=0.7\textwidth]{pic/fig-8-4-5.png}
% \caption{ }
\end{figure}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.5. 延迟微分方程 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：延迟微分方程经常出现在哪些学科中？

\item  答：控制论、生物数学等。

\item  问：有关延迟微分方程的文献有哪些？

\begin{itemize}
\item  Erneux, Thomas, Applied Delay Differential Equations, Springer 2009. 
\item  Driver, R.D., Ordinary and Delay Differential Equations, Springer 1997. 
\end{itemize}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.5.1. 模型方程 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：什么是延迟微分方程？
\item  答：未知函数在每个时刻的导数值由这个函数在过去时刻的值确定。初始值是零时刻之前的一段时间内的函数值。
\begin{eqnarray*}
\left\{\begin{array}{rcl}
\frac{dx}{dt}(t) &=& x(t-\tau), \,\, t>0, \\ 
x(t) &=& x_0(t), \,\, t\in [-\tau,0]. 
\end{array}\right. 
\end{eqnarray*}

\item  问：给定上述方程，解释在原点的单侧导数的不同含义，
\begin{eqnarray*}
\frac{dx}{dt}\Big{|}_{t\to 0+} &=& x(-\tau), \\
\frac{dx}{dt}\Big{|}_{t\to 0-} &=& \frac{dx_0}{dt}\Big{|}_{t\to 0-}. 
\end{eqnarray*}


\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.5.1. 一个教学例子  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：设 $x(t)$ 是定义在 $[-1,\infty)$ 的函数，满足下述条件，求解这个方程， 
\begin{eqnarray*}
\left\{\begin{array}{rcl}
\frac{dx}{dt}(t) &=& x(t-1), \,\, t>0, \\ 
x(t) &=& 1, \,\, t\in [-1,0]. 
\end{array}\right. 
\end{eqnarray*}

\item  答：

\begin{eqnarray*}
x(t) =
\left\{\begin{array}{ll}
1, & -1\le t\le 0, \\ 
1+t, & 0< t\le 1, \\ 
(3+t^2)/2, & 1< t\le 2, \\
\cdots, & \cdots \\ 
\end{array}\right. 
\end{eqnarray*}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.5.2. 更一般的方程及其数值解 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：有关延迟微分方程的数值解的文献？
\begin{itemize}
\item  Bellen, A. and Zennaro, M., Numerical Methods for Delay Differential Equations, Oxford, 2003.  
%\item  Bogacki, P. and Shampine, L. F., A3(2) pair of Runge-Kutta formulas, Applied Mathematics Letters,  2, 321-325. 
%\item  Shampine, L. F. and Thompson, S., Solving DDEs in Matlab, Appl. Num. Math. 37, 4, 441, 2001.
\item  Mackey, M. C. and Glass, L.. Pathological physiological conditions resulting from instabilities in physiological control system. Science, 197(4300):287-289. 1977.
%\item  Lang, R. and Kobayashi, K. , External optical feedback effects on semiconductor injection laser properties, IEEE J. Quantum Electron. 16, 347, 1980. 

\item  Flunkert, V. and Sch{\"o}ll, E., pydelay -- a python tool for solving delay differential equations, Arxiv, 2009.

\item  
\url{http://www.scholarpedia.org/article/Mackey-Glass_equation}

\end{itemize}


\end{itemize}

\end{frame}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.5.3. 逻辑斯谛延迟微分方程  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  逻辑斯谛微分方程：
\begin{eqnarray*}
\frac{dx(t)}{dt} = x(t)(1-x(t)). 
\end{eqnarray*}

\item  逻辑斯谛延迟微分方程：
\begin{eqnarray*}
\frac{dx(t)}{dt} &=& x(t)(1-x(t-\tau)),\,\, t\ge 0. \\
x(t) &=& x_0,\,\, -\tau\le t\le 0.
\end{eqnarray*}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.5.3. 种群动力学模型 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  种群动力学模型的一些文献：

\begin{itemize}
\item  Murray, J.D., Mathematical biology I. An Introduction, Springer, 2002.
\item  Erneux, Thomas, Applied Delay Differential Equations, Springer 2009. 
%\item  
%\item  
\end{itemize}



\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.5.3.   }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：对于 $\tau$ 的不同取值，逻辑斯谛延迟微分方程的解有什么不同？
\begin{itemize}
\item  当 $\tau=0$ 时
\item  当 $0 < \tau < \frac{1}{e}$ 时
\item  当 $\frac{1}{e} < \tau < \frac{\pi}{2}$ 时
\item  当 $\tau > \frac{\pi}{2}$ 时
\end{itemize}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.5.4. 麦克-格拉斯方程 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：什么是麦克-格拉斯方程？
\item  答：这是生物数学里的一个微分方程模型，
\begin{eqnarray*}
\left\{\begin{array}{rcl}
\frac{dx}{dt} (t) &=& a\frac{x(t-\tau)}{1+x(t-\tau)^m} - bx(t),\,\, t>0, \\ 
x(t) &=& x_0(t), \,\, t\in [-\tau,0]. 
\end{array}\right. 
\end{eqnarray*}

\item  问：有关延迟微分方程的数值解的软件？

\item  答：\url{http://pydelay.sourceforge.net}



\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.5.4. Introduction of pydelay }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  {\color{blue}pydelay} is a program which translates a system of delay differential equations (DDEs) into simulation C-code and compiles and runs the code (using {\color{blue}scipy weave}). This way it is easy to quickly implement a system of DDEs but you still have the speed of C. It is largely inspired by {\color{blue}PyDSTool}.

\item  The algorithm used is based on the Bogacki-Shampine method which is also implemented in Matlab's {\color{blue}dde23}.

{\color{blue}pydelay} is licensed under the MIT License.


\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.5.4.  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  代码1/3:

\begin{python}
import numpy as np
import pylab as pl
from pydelay import dde23

eqns = { 
    'x' : '0.25 * x(t-tau) / (1.0 + pow(x(t-tau),p)) -0.1*x' 
    }

params = {
    'tau': 15,
    'p'  : 10
    }
\end{python}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.5.4.  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  代码2/3:

\begin{python}
dde = dde23(eqns=eqns, params=params)

dde.set_sim_params(tfinal=1000, dtmax=1.0)

histfunc = {
    'x': lambda t: 0.5 
    } 

dde.hist_from_funcs(histfunc, 51)

dde.run()

\end{python}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.5.4.  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  代码3/3:

\begin{python}
sol1 = dde.sample(515, 1000, 0.1)
x1 = sol1['x']

sol2 = dde.sample(500, 1000-15, 0.1)
x2 = sol2['x']

pl.plot(x1, x2)
pl.xlabel('$x(t)$')
pl.ylabel('$x(t - 15)$')
\end{python}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.5.4. 麦克-格拉斯方程的一个周期解 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{figure}
\centering
\includegraphics[height=0.7\textheight, width=0.7\textwidth]{pic/mackey-glass.png}
% \caption{ }
\end{figure}


\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6. 随机微分方程 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：在生物数学和金融数学等领域，经常无法获得精确的观测值。找出一些研究随机演化规律的参考文献。

\begin{itemize}
\item  Gardiner, C. W., Handbook of Stochastic Methods, Fourth Edition, Springer, 2009. 
\item  Oksendal, B., Stochastic Differential Equations, Sixth Edition, Springer, 2003. 
\item  Kloeden, P. E. and Platen, E., Numerical Solution of Stochastic Differential Equations, Springer, 1992.
\item  Higham, D. J., An algorithmic introduction to numerical solution of stochastic differential equations, SIAM Review 43, 525-546. 
\end{itemize}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.1. 维纳过程  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：什么是维纳过程？（标准布朗运动）

\item  答：一系列连续依赖于 $t\in [0,T]$ 的随机变量 $W(t)$, 满足准则：
\begin{enumerate}
\item  $W(0)=0$ 的概率为 1.
\item  设 $0 \le s < t \le T$, 则增量 $W(t)-W(s)$ 是正态分布的随机变量，均值为零，方差为 $t-s$. 
\item  设 $0 \le s < t < u < v \le T$, 则增量 $W(t)-W(s)$ 与 $W(v)-W(u)$ 是相互独立的随机变量。
\end{enumerate}


\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.1.  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：构造一个维纳过程的样本路径。

\begin{python}
import numpy as np
import numpy.random as npr

T=1; N=500
t,dt=np.linspace(0,T,N+1,retstep=True)
dW=npr.normal(0.0,np.sqrt(dt),N+1)
dW[0]=0.0
W=np.cumsum(dW)

import matplotlib.pyplot as plt 
plt.plot(t,W)
plt.xlabel('t'); plt.ylabel('W(t)')
plt.title('Sample Wiener Process',weight='bold',size=12)
\end{python}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.1.  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{figure}
\centering
\includegraphics[height=0.7\textheight, width=0.7\textwidth]{pic/fig-8-6-1.png}
% \caption{ }
\end{figure}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.2. It{\^ o} 微积分 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：设维纳过程的增量 $dW(t) = W(t+\Delta t) - W(t)$. 求 $dW(t)$ 的均值和方差。

\item  答：因为 $dW(t) \sim N(0,\Delta t)$, 所以它的均值为零，方差为 $\Delta t$, 即
\begin{eqnarray*}
\mathbb{E}(dW) = 0, \,\, \mathbb{E}[(dW)^2] =dt.
\end{eqnarray*}


\item  注：设 $X$ 是随机变量，有时记 $\langle X \rangle$ 为其一阶矩即 $\mathbb{E}(X)$, 记 $\langle X^2 \rangle$ 为其二阶矩即 $\mathbb{E}(X^2)$.  

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.2. }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：什么是确定微分方程？答：设 $x(t)$ 是可微函数，
\begin{eqnarray*}
dx(t) &=& a(x(t)dt, \,\, x(0)=x_0,\\ 
x(t) - x_0 &=& \int_0^t a(x(s))ds. 
\end{eqnarray*}

\item  问：什么是随机微分方程？答：设 $X(t)$ 是随机过程，
\begin{eqnarray*}
dX(t) &=& a(X(t))dt + b(X(t))dW(t), \,\, X(0)=X_0, \\ 
X(t) - X_0 &=& \int_0^t a(X(s))ds + \int_0^t b(X(s))dW(s).  
\end{eqnarray*}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.2. It{\^o}公式 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：设 $X(t)$ 是随机过程， 满足下述随机微分方程，
\begin{eqnarray*}
dX(t) = a(X(t))dt + b(X(t))dW(t).  
\end{eqnarray*}
设 $Y(t)=f(X(t))$, 其中 $f$ 是可微函数，求 $Y(t)$ 满足的随机微分方程。

\item  答：按照泰勒展开，根据 $(dW)^2=dt$ 的原则，可得
\begin{eqnarray*}
dY &=& f\,'(X)dX + \frac{1}{2}f\,''(X)(dX)^2 + \cdots \\ 
&=& f\,'(X)[adt+bdW] + \frac{1}{2}f\,''(X) [adt+bdW]^2 + \cdots \\ 
&=& \left[ af\,'(X) + \frac{1}{2}b^2f\,''(X) \right] dt + bf\,'(X)dW + o(dt)
\end{eqnarray*}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.2. It{\^o}公式的例子  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：设 $X(t)$ 是维纳过程，设 $f(x)=e^{-x/2}$. 设 $Y=f(X)$. 
\begin{enumerate}
\item  求 $Y$ 满足的随机微分方程。
\item  求 $Y$ 的均值函数 $\langle Y \rangle$ 的解析解与数值解。
\end{enumerate}

\item  答：
\begin{enumerate}
\item  
在 It{\^o} 公式中，代入 $a=0,b=1$, 可得
\begin{eqnarray*}
dY &=& -\frac{1}{2}YdW + \frac{1}{2}\left(-\frac{1}{2}\right)^2Ydt, \\
Y(0) &=& \exp(-X(0)/2) = \exp(0)=1.
\end{eqnarray*}
\end{enumerate}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.2. 例子的解析解 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  答：
\begin{enumerate}
\item[2.]  均值函数的定义为 $\langle Y \rangle (t) := \mathbb{E}(Y(t))$. %求 $\langle Y \rangle (t) $ 的表达式。
在上式两边求期望，可得
\begin{eqnarray*}
d\langle Y\rangle = \frac{1}{8}\langle Y \rangle dt, \,\, \langle Y\rangle (0)=1.
\end{eqnarray*}
积分可得解析解
\begin{eqnarray*}
\langle Y \rangle (t) = \exp(t/8). 
\end{eqnarray*}
\end{enumerate}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.2. 例子的模拟解  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：模拟5000个维纳过程 $\{X(t),\,\, 0\le t\le 1\}$ 的样本路径，以此估计随机过程 $Y(t)=\exp(-X(t)/2)$ 的均值函数。

\item  代码1/2: 

\begin{python}
import numpy as np
import numpy.random as npr
T=1; N=1000; M=5000
t,dt=np.linspace(0,T,N+1,retstep=True)
dW=npr.normal(0.0,np.sqrt(dt),(M,N+1))
dW[ : ,0]=0.0
W=np.cumsum(dW,axis=1)
U=np.exp(- 0.5*W)
Umean=np.mean(U,axis=0)
Uexact=np.exp(t/8)
\end{python}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.2.  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  代码2/2: 
{\footnotesize 
\begin{python}
import matplotlib.pyplot as plt 
fig=plt.figure(); ax=fig.add_subplot(111)
ax.plot(t,Umean,'b-',label='mean of %d paths' % M)
ax.plot(t,Uexact,'r-.', label='exact mean function')
for i in range(5):
    ax.plot(t,U[i,:],'g--')

ax.set_xlabel('t')
ax.set_ylabel('U')
ax.legend(loc='best')

maxerr=np.max(np.abs(Umean-Uexact))
print('With %d paths and %d intervals the max error is \  
        %g' % (M,N,maxerr))
\end{python}
}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.2. 解析解与模拟解的图像  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{figure}
\centering
\includegraphics[height=0.7\textheight, width=0.7\textwidth]{pic/fig-8-6-2-b.png}
% \caption{ }
\end{figure}

\end{frame}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.3. Riemann-Stieltjes 积分  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：什么是黎曼-斯蒂尔杰斯积分？

\item  答：将积分区间 $[a,b]$ 分成一些小区间 $a=t_0<t_1<\cdots<t_N=b$, 在每个小区间上任取 $\tau_k\in [t_k,t_{k+1}]$, 如果下述和式当区间最大宽度趋于零时，存在与分划和选取点无关的极限，则称其为 RS 积分，
$$ \int_{a}^{b} f(t) dg(t) = \lim\limits_{N\to\infty} \sum\limits_{k=0}^{N-1} f(\tau_k) [g(t_{k+1})-g(t_k)]. $$


\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.3.  It{\^o} 积分 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：什么是伊藤积分？

\item  答：在RS积分中，取 $g(t)=W(t)$ 为布朗运动，然后在每个小区间，取左端点 $\tau_k=t_k$, 这样得到的RS和式是一个随机变量，极限理解为这个随机变量序列在二阶矩收敛意义下的极限，
$$ \int_{T_1}^{T_2} f(t) dW(t) = \lim\limits_{n\to\infty} \sum\limits_{k=0}^{n-1} f(\tau_k) [W(t_{k+1})-W(t_k)]. $$

\item  注：称随机变量序列 $\{\xi_n,\,\, n\ge 1\}$ 二阶矩收敛于随机变量 $\xi$, 是指
$$\lim\limits_{n\to\infty} \mathbb{E}(\xi_n) =\mathbb{E}(\xi), \,\,\, \lim\limits_{n\to\infty} \text{Var}(\xi_n-\xi) =0. $$


\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.3. Stratonovich 积分 }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：什么是斯特拉托诺维奇积分？

\item  答：在RS积分中，取 $g(t)=W(t)$ 为布朗运动，然后在每个小区间，取中间点 $\tau_k=(t_k+t_{k+1})/2$, 
$$ \int_{T_1}^{T_2} f(t) dW(t) = \lim\limits_{n\to\infty} \sum\limits_{k=0}^{n-1} f(\tau_k) [W(t_{k+1})-W(t_k)]. $$


\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.3. 一个例子  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：固定参数 $\lambda\in [0,1]$, 在每个小区间 $[t_k,t_{k+1}]$ 取固定比例的点 
$$\tau_k = (1-\lambda)t_k + \lambda t_{k+1}, $$ 
\begin{enumerate}
\item  计算 RS 和式的均值与方差，验证其二阶矩极限是下述右边，
$$\int_{T_1}^{T_2} W(t)dW(t) = \frac{1}{2} \left[ W(T_2)^2 - W(T_1)^2 \right] + \left( \lambda-\frac{1}{2} \right)(T_2-T_1). $$
\item  使用模拟样本路径，验证上述等式。
\end{enumerate}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.3.  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  答：
\begin{enumerate}
\item  记 $W_k := W(t_k)$, $W_{k+1} := W(t_{k+1})$, 以及 $W_{\tau} := W(\tau_{k})$ (not a good notation!), 

\item  则 RS 和式的每一项为 
\begin{eqnarray*}
 W_{\tau} (W_{k+1} - W_k) 
&=& \frac{1}{2} (W_{k+1}^2 - W_k^2) - \frac{1}{2} (W_{k+1} - W_k)^2 \\
&& + (W_\tau - W_k)^2 + (W_{k+1} - W_\tau)(W_\tau - W_k)
\end{eqnarray*}

\item  求和，根据 $(dW(t))^2=dt$ 可得所证等式的右边。

\item  将上一节的代码略加修改，即可验证。
\end{enumerate}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.4. 随机微分方程的数值求解  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：求下述随机微分方程的数值解， 
\begin{eqnarray*}
dX(t) = a(X(t))dt + b(X(t))dW(t), \,\, X(0)=X_0. 
\end{eqnarray*}

\item  答：将时间区间 $[0,T]$ 等分成 $N$ 个小区间，即 $t_k=\frac{kT}{N}$. 记 $X_{k}=X(t_k)$, 
\begin{eqnarray*}
X_{k+1} = X_k + \int_{t_k}^{t_{k+1}} a(X(s))ds + \int_{t_k}^{t_{k+1}} b(X(s))dW(s).  
\end{eqnarray*}
问题转化为如何近似计算这两个积分。
\begin{enumerate}
\item  欧拉-丸山数值解法
\item  米尔斯坦数值解法
\end{enumerate}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.4.  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：什么是 Euler-Maruyama method 和 Milstein method?
\item  答：是将 It{\^o} 积分用下述来近似的方法，
\begin{eqnarray*}
X_{k+1} &=& X_k + a(X_k)\Delta t + b(X_k)(W_{k+1}-W_k), \\  
X_{k+1} &=& X_k + a(X_k)\Delta t + b(X_k)(W_{k+1}-W_k) \\
&&\,\,\,\,\,\,\,\, + \frac{1}{2}b(X_k)b'(X_k)[(W_{k+1}-W_k)^2-\Delta t].
\end{eqnarray*}


\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.4.  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：欧拉-丸山方法和米尔斯坦方法的原理是什么？

\item  答：这两个数值方法都是为了估计下述积分
\begin{eqnarray*}
X_{k+1} = X_k + \int_{t_k}^{t_{k+1}} a(X(s))ds + \int_{t_k}^{t_{k+1}} b(X(s))dW(s).  
\end{eqnarray*}
%设 $a(X(t))$ 和 $b(X(t))$ 
考虑一般的 $Y(X(t))$. 将其微分展开，可得
\begin{eqnarray*}
dY &=& Y'(X)dX + \frac{1}{2} Y''(X)^2(dX)^2 \\
&=& Y'(X) (a(X)dt+b(X)dW) + \frac{1}{2}Y''(X)^2 (a(X)dt+b(X)dW))^2 \\
&=& \mathcal{L}[Y]dt +\mathcal{M}[Y]dW + \text{高阶项},
\end{eqnarray*}
其中的记号为
$%\begin{eqnarray*}
\mathcal{L}[Y] = a(X)Y'(X) + \frac{1}{2}Y''(X)^2b(X)^2,  \,\,\,\,
\mathcal{M}[Y] = b(X)Y'(X).
$%\end{eqnarray*}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.4.  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  答：将这个一般的 $Y(X(t))$ 的结论用于 $a(X(t))$ 和 $b(X(t))$, 可得
\begin{eqnarray*}
a(X(s)) &=& a(X(t_k)) + \int_{t_k}^s \mathcal{L}[a]dt + \int_{t_k}^s \mathcal{M}[a] dW, \\
b(X(s)) &=& b(X(t_k)) + \int_{t_k}^s \mathcal{L}[b]dt + \int_{t_k}^s \mathcal{M}[b] dW. 
\end{eqnarray*}
代入所需估计的两个积分，得到
\begin{eqnarray*}
X_{k+1} &=& X_k + a(X_k)(t_{k+1}-t_k) + b(X_k)(W(t_{k+1}) - W(t_k)) \\
&& + \int_{t_k}^{t_{k+1}} \left[ \int_{t_k}^s b(X(\tau))b'(X(\tau)) dW(\tau) \right]  dW(s) + \text{高阶项}.
\end{eqnarray*}
保留前三项得欧拉-丸山方法，保留前四项得米尔斯坦方法。

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.4.  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：考虑 Black-Scholes 方程 
\begin{eqnarray*}
dX(t) = \lambda X(t)dt + \mu X(t)dW(t), \,\, X(0)=X_0. 
\end{eqnarray*}

\begin{enumerate}
\item  写出解析解。
\item  计算数值解。
\end{enumerate}

\item  答：解析解如下，这个随机过程称为几何布朗运动，
$$X(t) = X_0 \exp\left[ \left(\lambda - \frac{1}{2} \mu^2\right) t + \mu W(t) \right]. $$
对维纳过程的每条路径，得到解析解的路径和数值解的路径，下面比较的是这两条路径的区别。



\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.4.  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  代码1/4：

\begin{python}
import numpy as np
import numpy.random as npr

#set up grid
T=1.0
N=30
t,dt=np.linspace(0,T,N+1,retstep=True)

#Get Brownian motion
dW=npr.normal(0.0,np.sqrt(dt),N+1)
dW[0]=0.0
W=np.cumsum(dW)
\end{python}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.4.  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  代码2/4：

\begin{python}
#Equation parameters and functions
lamda=2.0
mu=1.0
Xzero=1.0
def a(X): return lamda*X
def b(X): return mu*X
def bd(X): return mu*np.ones_like(X)

#Analytic solution
Xanal=Xzero*np.exp((lamda-0.5*mu*mu)*t+mu*W)
\end{python}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.4.  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  代码3/4：

\begin{python}
#Milstein solution
Xmil=np.empty_like(t)
Xmil[0]=Xzero
for n in range(N):
    Xmil[n+1]=Xmil[n] + dt*a(Xmil[n]) + dW[n+1]*b(Xmil[n]) \
    +0.5*(b(Xmil[n])*bd(Xmil[n])*(dW[n+1]**2-dt))
\end{python}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.4.  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  代码4/4：

\begin{python}
import matplotlib.pyplot as plt
plt.ion()
plt.plot(t,Xanal,'b-',label='analytic')
plt.plot(t,Xmil,'g-.',label='Milstein')
plt.legend(loc='best')
plt.xlabel('t')
plt.ylabel('X(t)')
\end{python}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.4. 解析解与数值解的比较  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{figure}
\centering
\includegraphics[height=0.7\textheight, width=0.7\textwidth]{pic/fig-8-6-4-milstein.png}
% \caption{ }
\end{figure}

\end{frame}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.4.  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：什么是强收敛的阶？
\item  答：将数值解 $\{X_0,X_1,\cdots,X_n,\cdots,X_N\}$ 与解析解 $\{X(\tau),t_0\le \tau \le t_N\}$ 进行比较，若下述成立，则称这个方法有强收敛的阶数 $\gamma$, 
$$\mathbb{E} |X_n - X(\tau) | = O( (\Delta t)^{\gamma} ), \,\,\, \text{当} \Delta t\to 0. $$

\item  问：验证欧拉-丸山方法和米尔斯坦方法的强收敛的阶分别是 $1/2$ 和 $1$.
\item  答：
\begin{enumerate}
\item  理论证明。
\item  模拟验证。
\end{enumerate}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.4.  }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：什么是弱收敛的阶？
\item  答：将数值解 $\{X_0,X_1,\cdots,X_n,\cdots,X_N\}$ 与解析解 $\{X(\tau),t_0\le \tau \le t_N\}$ 进行比较，若下述成立，则称这个方法有弱收敛的阶数 $\gamma$, 
$$ | \mathbb{E} (X_n) - \mathbb{E} (X(\tau) | = O( (\Delta t)^{\gamma} ), \,\,\, \text{当} \Delta t\to 0. $$

\item  问：验证欧拉-丸山方法和米尔斯坦方法的弱收敛的阶都是1.
\item  答：
\begin{enumerate}
\item  理论证明。
\item  模拟验证。
\end{enumerate}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.4.   }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  问：模拟样本路径，研究欧拉-丸山方法和米尔斯坦方法的强收敛的阶。

\item  答：按定义，如果下述成立，则称数值解 $\{X_n\}$ 的强收敛的阶为 $\gamma$. 
\begin{eqnarray*}
\mathbb{E} |X_n - X(\tau) | = O( (\Delta t)^{\gamma} ), \,\,\, \text{当} \Delta t\to 0.
\end{eqnarray*}

画图验证这两个阶分别为 $\frac{1}{2}$ 与 $1$. 为直观呈现，取对数可得
\begin{eqnarray*}
\log \mathbb{E} |X_n - X(\tau) | \sim \gamma \log (\Delta t).
\end{eqnarray*}
因此 $\gamma$ 是横坐标为 $\log(\Delta t)$, 纵坐标为 $\log(Error)$ 的图像的斜率。


\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.4.   }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  代码1/6：

\begin{python}
import numpy as np
import numpy.random as npr

#Problem definition
M=1000  #number of paths sampled
P=6  #number of discretizations
T=1  #endpoint of time interval
N=2**12  #finest grid size
dt=1.0*T/N

#Problem parameters
lamda=2.0
mu=1.0
Xzero=1.0
\end{python}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.4.   }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  代码2/6：

\begin{python}
def a(X): return lamda*X
def b(X): return mu*X
def bd(X): return mu*np.ones_like(X)

#build the Brownian paths
dW=npr.normal(0.0, np.sqrt(dt),(M,N+1))
dW[:,0]=0.0
W=np.cumsum(dW,axis=1)

#build the exact solutions at the end of the paths
ones=np.ones(M)
Xexact=Xzero*np.exp((lamda-0.5*mu*mu)*ones+mu*W[:,-1])
Xemerr=np.empty((M,P))
Xmilerr=np.empty((M,P))
\end{python}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.4.   }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  代码3/6：

\begin{python}
#loop over refinements
for p in range(P):
    R=2**p
#    L=int(N/R)
    L=2**(12-p);Dt=R*dt
    Xem=Xzero*ones;Xmil=Xzero*ones
    Wc=W[:,::R]
    for j in range(L):
        deltaW=Wc[:,j+1]-Wc[:,j]
        Xem+=Dt*a(Xem)+deltaW*b(Xem)
        Xmil+=Dt*a(Xmil)+deltaW*b(Xmil)+0.5*b(Xmil)*bd(Xmil)*(deltaW**2-Dt)
    Xemerr[:,p]=np.abs(Xem-Xexact)
    Xmilerr[:,p]=np.abs(Xmil-Xexact)
\end{python}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.4.   }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  代码4/6：

\begin{python}
#do some plotting
import matplotlib.pyplot as plt
Dtvals=dt*np.array([2**p for p in range(P)])
lDtvals=np.log10(Dtvals)
Xemerrmean=np.mean(Xemerr,axis=0)
plt.plot(lDtvals,np.log10(Xemerrmean),'bo')
plt.plot(lDtvals,np.log10(Xemerrmean),'b:',label='EM actual')
plt.plot(lDtvals,0.5*np.log10(Dtvals),'b-.',label='EM theoretical') 
\end{python}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.4.   }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  代码5/6：

\begin{python}
Xmilerrmean=np.mean(Xmilerr,axis=0)
plt.plot(lDtvals,np.log10(Xmilerrmean),'go')
plt.plot(lDtvals,np.log10(Xmilerrmean),'g:',label='Mil actual')
plt.plot(lDtvals,np.log10(Dtvals),'g-.',label='Mil theoretical') 
plt.legend(loc='best')
#plt.xlabel('log of Delta t',size=12)
#plt.ylabel('log of expectation of difference',size=12)
plt.xlabel(r'$\log_{10}\Delta t$',size=16) 
plt.ylabel(r'$\log_{10}\left(\langle|X_n-X(\tau)|\rangle\right)$',size=16)
\end{python}

\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.4.   }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{itemize}

\item  代码6/6：

\begin{python}
emslope=((np.log10(Xemerrmean[-1]) - np.log10(Xemerrmean[0]))/
         (lDtvals[-1]-lDtvals[0]))
print('Emprical EM slope is %g' %emslope)
milslope=((np.log10(Xmilerrmean[-1]) - np.log10(Xmilerrmean[0]))/
          (lDtvals[-1]-lDtvals[0]))
print('Emprical mil slope is %g' %milslope)
\end{python}


\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile=singleslide]{8.6.4.   }
%\begin{frame}{}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{figure}
\centering
\includegraphics[height=0.7\textheight, width=0.7\textwidth]{pic/fig-8-6-4-strong-convergence.png}
% \caption{ }
\end{figure}

\end{frame}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\begin{frame}[fragile=singleslide]{1.20. }
\begin{frame}{参考文献}

\vspace{-0.4cm}\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%每页详细内容

\begin{thebibliography}{99}

\bibitem{stewart-en} John M. Stewart. \emph{Python for Scientists}. Second Edition. Cambridge University Press. 2017. 
\bibitem{stewart-cn} 约翰.M.斯图尔特(著). 江红.余青松(译). \emph{Python科学计算}，机械工业出版社，2019年8月第1版。

\bibitem{sauer-en} Timothy Sauer. \emph{Numerical Analysis}. Third Edition. Pearson. October 2017. 
\bibitem{sauer} Timothy Sauer(著).裴玉茹.马赓宇(译). \emph{数值分析}. 机械工业出版社. 2018年8月第1版.     


\end{thebibliography}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\end{document}


